# include "scene.h"

#ifndef CALL_SOLVE_POISSON_OMP

# include "discr.h"
# include <stdio.h>
# include <malloc.h>
# include <math.h>
/*
# include <time.h>
# include <sys/timeb.h>
# include <omp.h>*/

extern void press_bc(const struct grid *g, short tflag, unsigned short iRK);

extern long int solve_poisson_old(struct grid *g, const double ***d, unsigned short iRK, double ***fspace)
  {
   int cnt;
   short i, j, k, nr, nth, nz, im, jm, km, i1, j1, k1, i2, j2, k2, ii, jj, kk, stop_cnt;
   double dpm, dp, ***p, *r, *th, *z, prf, ***pmemloc;
   char flag;

// struct timeb t1,t2;

#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'solve_poisson()'...\n");
#  endif

   r = g->r;   z = g->z;   th = g->th;   nr = g->nr;   nth = g->nth;   nz = g->nz;

   prf = PRF_INIT; stop_cnt = 0;

   pmemloc = g->p;  p = g->p;
   flag = 'o';
   for(cnt = 1; cnt <= POISS_CNT_MAX; cnt++)
      {
#      ifdef PPE_ITTR_JACOBI
       if(flag == 'o'){ p = pmemloc; g->p = fspace; flag = 'e';}
       else{ p = fspace; g->p = pmemloc; flag = 'o';}
#      endif

       if(cnt%2 == 1) { k1 = NzG; k2 = NzG+nz; kk = 1; }
       else { k1 = NzG+nz-1; k2 = NzG-1; kk = -1; }

       for(dpm = 0.0, im = jm = km = -1, k = k1; k != k2; k += kk)
          {
           if(cnt%2 == 1) { i1 = NrG; i2 = NrG+nr; ii = 1; }
           else { i1 = NrG+nr-1; i2 = NrG-1; ii = -1; }

           for(i = i1; i != i2; i += ii)
               {
                if(cnt%2 == 1) { j1 = NthG; j2 = NthG+nth; jj = 1; }
                else { j1 = NthG+nth-1; j2 = NthG-1; jj = -1; }

                for(j = j1; j != j2; j += jj)
                   {
                    dp = P_ijk_NEW - p[i][j][k];

                    if((dp<0 ? -dp:dp) > (dpm<0 ? -dpm:dpm))
                         {
                          dpm = dp;
                          im = i; jm = j; km = k;
                         }
                    g->p[i][j][k] = p[i][j][k] + prf*dp;
                   }
               }
          }

       press_bc(g, 0, iRK);

       if( (cnt == 1) || (cnt%DISP_FREQ == 0) )
            {
             fprintf(g->dptr, "=> %d: dpm = %E(%d, %d, %d)\n", cnt, dpm, im, jm, km);
             fflush(g->dptr);
            }

       if( ((dpm < 0) ? -dpm:dpm) <= PRESS_DELTA) ++stop_cnt;
       else stop_cnt = 0;

       if(stop_cnt == POISS_CNT_MIN) break;
      }

/* If the new pressure field is located at "fspace", copy it to the memory located at _p. */
   if(flag == 'e')
      {
       for(i = 0; i < 2*NrG+nr; i++) for(j = 0; j < 2*NthG+nth; j++) for(k = 0; k < 2*NzG+nz; k++)
       pmemloc[i][j][k] = fspace[i][j][k];

       g->p = pmemloc;
      }

   flag = (stop_cnt == POISS_CNT_MIN) ? 'Y':'N';

   fprintf(g->dptr, "N_ittr = %d; dpm = %E(%d, %d, %d); Conv. = %c\n", cnt, dpm, im, jm, km, flag);
   fflush(g->dptr);

   fprintf(g->tptr, "N_ittr = %d; dpm = %E(%d, %d, %d); Conv. = %c; ", cnt, dpm, im, jm, km, flag);
   fflush(g->tptr);

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'solve_poisson()' ends.\n");
#  endif

   return cnt;
  }

#endif
